'''
Created on Mar 17, 2010

@author: oabalbin
'''

import numpy as np
#import scipy as sp
#from scipy import stats
#import matplotlib.pyplot as plt
#import matplotlib.mlab as mlab
# import rpy2 modules
import rpy2.robjects as robjects
import rpy2.robjects.numpy2ri
from rpy2.robjects.packages import importr


def load_common_Rlibraries():
    """
    This function load some common libraries in R that can be used for other of the functions
    this package
    """
    r = robjects.r
    r("library(robust)")


def calc_robust_correlation(matData,libloaded):
    """
    It calculates a Robust correlation for matData. 
    It uses R function for covRob. This function is especially
    important when the data have many outliers that can drive the correlation.
    """
    
    r = robjects.r
    if not libloaded:
        r("library(robust)")
    
    naexclude = r["na.exclude"]
    covRob = r["covRob"]
    corr_results =  covRob(matData,corr="TRUE", estim="pairwiseQC", na_action=naexclude)
    corrMat = np.array(corr_results[1])
    
    cor_val = corrMat[0,1]
    
    return cor_val


def calc_matrix_correlation(corrtype, myData, matrices=False, myDatay=[]):
    """
    It calculates the correlation between the columns of a matrix.
    The pearson correlation coefficient assumes that all observations where
    derived form the same normal distribution.  
    cor(x, y = NULL, use = "all.obs",
     method = c("pearson", "kendall", "spearman"))

    """
    parametric = {'Kend':'kendall','Pear':'pearson', 'Sper':'spearman'}[corrtype]
    
    r = robjects.r
    rcorr = r["cor"]
    # calc correlation matrix
    if matrices:
        corr_matrix = np.array( rcorr(myData, myDatay, method=parametric, use="pairwise") ) # use="pairwise", use="complete"
    else:
        corr_matrix = np.array( rcorr(myData,method=parametric, use="pairwise") ) # use=
    
    return corr_matrix


def plot_correlation_matrix(corrMat, pltname,tablefile): #myrownames, mycolnames
    """
    Plot a correlation matrix using corrgram 
    """
    r = robjects.r
    r("library(corrgram)")
    corrgram = r["corrgram"]
    
    panel_shade = r["panel.shade"]
    panel_pts = r["panel.pts"]
    panel_pie = r["panel.pie"]
    panel_txt = r["panel.txt"]
    dataframe = r["data.frame"]
    readtable = r["read.table"]
    panel_minmax=r["panel.minmax"]
    
    corrMat2 = readtable(tablefile,header=True, row_names=1,sep="\t")
    
    #corrMat2 = dataframe(corrMat,row_names=1, col_names=mycolnames)
    
    #pltname = '/exds/users/oabalbin/projects/rnaseq_data/noveltu_expression_data/correlation_plot.png'
    
    grdevices = importr('grDevices')
    grdevices.png(file=pltname, width=15000, height=10000)
    r.png(pltname, width=15000, height=10000)
    
    corrgram(corrMat2, order=False, lower_panel=panel_shade, upper_panel=panel_pie, text_panel=panel_txt,
      main="EZH2 targets vs TUs in PC2/PC1 Order") 
    
    grdevices.dev_off()



def plot_correlation_matrix2(corrMat, pltname, tablefile):
    """
    """
    r = robjects.r
    r("source('/home/oabalbin/workspace/Rcode/Rplots.R')")
    
    dataframe = r["data.frame"]
    readtable = r["read.table"]
    asmatrix=r["as.matrix"]
    
    #corrMat2 = dataframe(corrMat,row_names=myrownames, col_names = mycolnames)
    corrMat2 = readtable(tablefile,header=True, row_names=1,sep="\t")
    corrMat2 = asmatrix(corrMat2)
    r.plot_correlation_matrix1(corrMat2, pltname)


def calc_hypergeometric_test(drawn_white_balls,total_white_balls, total_balls_drawn,total_black_balls, FDR):
    """
    It calculates an hypergeometric test using the parameters given
    phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
    Arguments
    x, q     vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
    m     the number of white balls in the urn.
    n     the number of black balls in the urn.
    k     the number of balls drawn from the urn.
    log, log.p     logical; if TRUE, probabilities p are given as log(p).
    lower.tail     logical; if TRUE (default), probabilities a    
    """
    
    r = robjects.r
    r("library(stats)")
    r("library(qvalue)")
    
    rphyper=r["phyper"]
    rqvalue=r["qvalue"]
    #rlambda=r["lambda"] 
    
    prob_enriched = rphyper(drawn_white_balls - 1, total_white_balls, total_black_balls, total_balls_drawn, lower_tail=False);
    
    #myqvalue = rqvalue(prob_enriched,fdr_level=FDR)
    #print np.array(myqvalue)
    
    return np.array(prob_enriched)[0] #, np.array(myqvalue)[0]#np.array(prob_enriched)[0], np.array(myqvalue)[0]


########## This deprecated use R script instead.
##########
def apply_randomForest(mF,samples,classes, Rcodefolder):
    r = robjects.r
    

    r("setwd('%s')" % Rcodefolder)
    
    r("library(randomForest)")
    r("source('mydataFrame.txt')")
    r("source('myRandomForest.txt')")
    r("source('tuneRandomForest.txt')")
    "mydataFrame"       "runRandomForest"   "tuneRandomForest2"

    mydataFrame = r['mydataFrame']
    runRandomForest = r["runRandomForest"]
    tuneRandomForest = r["tuneRandomForest2"]
    
    # Create a data frame of the mF (should be the transpose of mF)
    # and sample and ETS information
    myData = mydataFrame(mF,samples,classes)
    myData_rf = tuneRandomForest()
    
    r("tuneRandomForest2('%f','%d')" % myData,1000)


'''
drawn_white_balls,total_white_balls, total_balls_drawn,total_black_balls=4,5,10,45
calc_hypergeometric_test(drawn_white_balls,total_white_balls, total_balls_drawn,total_black_balls) 


mymat = np.loadtxt('/home/oabalbin/projects/networks/output/2010_03_20_21_13/drugsen/nci60/fact_ERG_vectors.txt',dtype=float, delimiter='\t')
mymat = np.transpose(mymat)
corval = calc_robust_correlation(mymat)
print corval
'''


'''
myData = np.loadtxt('/home/oabalbin/projects/networks/output/2010_03_20_21_13/bfrm/mF.txt',dtype=float, delimiter='\t')

#samples = np.loadtxt('/home/oabalbin/projects/networks/output/2010_03_15_00_41/bfrm/2010_03_15_00_41_dataset.sample.txt',dtype=str, delimiter='\t')
#classes = np.loadtxt('/home/oabalbin/workspace/signatures/trunk/signatures/stats/Rcode/ETS_class.txt',dtype=str, delimiter='\t')
#print samples, classes#
#Rcodefolder='/home/oabalbin/workspace/signatures/trunk/signatures/stats/Rcode/'
myDatar = np.random.uniform(0,1,(myData.shape[0],(myData.shape[1]-5) ) )

mat2mask = np.ma.make_mask_none((myDatar.shape[0],myDatar.shape[1]))    
mat2_raw = np.ma.array(myDatar,mask=mat2mask)
mat2_raw.mask[5:6,2:4] = True
#mat2_raw = np.ma.masked_where(myDatar < 0.01 , mat2_raw, copy=True)
print mat2_raw
#myDatar[5:7,2:4] = np.nan
print myDatar
print myData.shape
print myDatar.shape
corrtype="Sper"
var = range(myDatar.shape[1])
var = map(str, var)
var = np.array(var)
labels = range(myData.shape[1])
labels = map(str, labels)
labels = np.array(labels)

pltname = '/exds/users/oabalbin/projects/rnaseq_data/noveltu_expression_data/correlation_plot.png'

correlation_matrix = calc_matrix_correlation(corrtype, myData,True, mat2_raw.filled(fill_value=np.nan))

print correlation_matrix
print correlation_matrix.shape
print var
print labels
plot_correlation_matrix(correlation_matrix, pltname, labels, var)
#apply_randomForest(mF,samples,classes, Rcodefolder)

'''    

